/// Script to produce relative risk graph and estimates using nlcom command and coefplot package
/// Fabien Cottier, March 2019
/// CLIM DATA


///////// Graph immediate clim effect only (CLIM Y0)
* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est_clim	

	local modelsT_L0 ""
	local modelsT_L1 ""
	local modelsT_L2 ""

	
	local modelsP_L0 ""
	local modelsP_L1 ""
	local modelsP_L2 ""
	
	local T=-0.75
	local P=-1.25
	
		
	* TEMP / Loop to create multiple nlcom statement
	local j1=3
	while `T'<=1 {

		if "`2'"=="linear"{
			local modelsT `"`modelsT' (_nl_`j1':_b[stdw10ma_tempy]*`T')"'
		}
		
		else if "`2'"=="quadratic"{
			local modelsT `"`modelsT' (_nl_`j1':_b[stdw10ma_tempy]*`T'+_b[stdw10ma_tempy#stdw10ma_tempy]*(`T')^2)"'
		}
		
		local T=`T'+.25
		local j1=`j1'+1
	}
		
	* PRECIP / Loop to create multiple nlcom statement
	local j2=1
	while `P'<=1.25  {

		if "`2'"=="linear"{
			local modelsP `"`modelsP' (_nl_`j2':_b[stdw10ma_precipy]*`P')"'
		}
		
		else if "`2'"=="quadratic"{
			local modelsP `"`modelsP' (_nl_`j2':_b[stdw10ma_precipy]*`P'+_b[stdw10ma_precipy#stdw10ma_precipy]*(`P')^2)"'
		}
		
		local P=`P'+.25
		local j2=`j2'+1
	}
	
	
	* estimates relative risk using nlcom
	est resto `1' 
	local df_est=e(df_r)
	nlcom `modelsT', post df(`df_est')
	estimates store est_rr_T

	est resto `1' 
	nlcom `modelsP', post
	estimates store est_rr_P
	
	* produce plot using coefplot
	coefplot (est_rr_P), bylabel( ) || ///
		(est_rr_T), bylabel( ) ||, ///
		eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+75%") ///
		ytitle(Relative change in irr. migration (%)) vertical recast(line) ciopts(recast(rarea)) yline(1)
		addplot 1: , b1title("Std Precipitation deviation") norescaling  xlabel(1 "-1.25" 2 "-1" 3 "-.75" 4 "-.5" 5 "-.25" 6 "0" 7 ".25" 8 ".5" 9 ".75" 10 "1" 11 "1.25")
		addplot 2: , b1title("Std Temperature deviation") norescaling xlabel(1 "-1.25" 2 "-1" 3 "-.75" 4 "-.5" 5 "-.25" 6 "0" 7 ".25" 8 ".5" 9 ".75" 10 "1" 11 "1.25")
end



///////// Split sample graph clim (clim Y0)


* write programe for relative risk estimation // margins cannot produce rr estimates for cont. variables
program define rr_est_clim_Sp	

	local modelsT ""
	local modelsP ""
	
	local T=-0.75
	local P=-1.25
	
	* only to be executed if temperature estimates called
	if "`3'"=="temp"{
	
		* TEMP / Loop to create multiple nlcom statement
		local j1=3
		while `T'<=1 {

			if "`4'"=="linear"{
				local modelsT_L0 `"`modelsT_L0' (_nl_`j1':_b[stdw10ma_tempy]*`T')"'
			}
			
			else if "`4'"=="quadratic"{
				local modelsT_L0 `"`modelsT_L0' (_nl_`j1':_b[stdw10ma_tempy]*`T'+_b[stdw10ma_tempy#stdw10ma_tempy]*(`T')^2)"'
			}
			
		local T=`T'+.25
		local j1=`j1'+1
		}
		
		* estimates relative risk using nlcom
		est resto `1' 
		local df_est=e(df_r)
		nlcom `modelsT_L0', post df(`df_est')
		estimates store est_rr_TL0H
		
		est resto `2' 
		local df_est2=e(df_r)
		nlcom `modelsT_L0', post df(`df_est2')
		estimates store est_rr_TL0L
		
		* produce plot using coefplot
		coefplot (est_rr_TL0H), bylabel(High Agriculture - Immediate effect) || ///
			(est_rr_TL0L), bylabel(Low Agriculture - Immediate effect) || , ///
			eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+75%") ///
			ytitle(Relative change in irr. migration (%)) vertical recast(line) ciopts(recast(rarea)) yline(1) byopts(row(3)) ///
			xtitle("Std Temperature deviation") xlabel(1 ".75" 2 "-.5" 3 "-.25" 4 "0" 5 ".25" 6 ".5" 7 ".75" 8 "1")
			
	}
	
	* only to be executed if precipitation estimates called
	
	if "`3'"=="precip"{
	
		* PRECIP / Loop to create multiple nlcom statement
		local j2=1
		while `P'<=1.25  {

			if "`4'"=="linear"{
				local modelsP_L0 `"`modelsP_L0' (_nl_`j2':_b[stdw10ma_precipy]*`P')"'
			}
			
			else if "`4'"=="quadratic"{
				local modelsP_L0 `"`modelsP_L0' (_nl_`j2':_b[stdw10ma_precipy]*`P'+_b[stdw10ma_precipy#stdw10ma_precipy]*(`P')^2)"'
			}
			
		local P=`P'+.25
		local j2=`j2'+1
		}
		
		*estimates relative risk using nlcom
		est resto `1'
		local df_est=e(df_r)
		nlcom `modelsP_L0', post df(`df_est')
		estimates store est_rr_PL0H

		est resto `2' 
		local df_est2=e(df_r)
		nlcom `modelsP_L0', post df(`df_est2')
		estimates store est_rr_PL0L

		* produce plot using coefplot
		coefplot (est_rr_PL0H), bylabel(High Agriculture - Immediate effect) || ///
			(est_rr_PL0L), bylabel(Low Agriculture - Immediate effect) || , ///
			eform yscale(range(.75 1.75)) ylabel(.75 "-25%" 1 "0" 1.25 "+25%" 1.5 "+ 50%" 1.75 "+75%") ///
			ytitle(Relative change in irr. migration (%)) vertical recast(line) ciopts(recast(rarea)) yline(1) byopts(row(3)) ///
			xtitle("Std Precipitation deviation") xlabel(1 "-1.25" 2 "-1" 3 "-.75" 4 "-.5" 5 "-.25" 6 "0" 7 ".25" 8 ".5" 9 ".75" 10 "1" 11 "1.25")

	}
	
	
end



